Entanglements in Quiescent and Sheared Polymer Melts 
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We visualize entanglements in polymer melts using molecular dynamics simulation. A bead at 
an entanglement interacts persistently for long times with the non-bonded beads (those excluding 
the adjacent ones in the same chain). The interaction energy of each bead with the non-bonded 
beads is averaged over a time interval r much longer than microscopic times but shorter than the 
onset time of tube constraints r c . Entanglements can then be detected as hot spots consisting of 
several beads with relatively large values of the time-averaged interaction energy. We next apply a 
shear flow with rate much faster than the entangle motion. With increasing strain the chains take 
zigzag shapes and a half of the hot spots become bent. The chains are first stretched as a network 
but, as the bends approach the chain ends, disentanglements subsequently occur, leading to stress 
overshoot observed experimentally. 



PACS numbers: 83.10.Rs, 83. 50. Ax, 61.25.Hq 



I. INTRODUCTION 

The dynamics of dense polymer melts has been a chal- 
lenging subject in current polymer physics [1-5]. While 
the near-equilibrium dynamics of short chains with N < 
N e can be reasonably well described by the simple Rouse 
model, the dynamics of very long chains with N > N 
has not yet been well understood on the microscopic level 
since it is governed by complicated entanglement effects. 
Here N is the polymerization index or the bead number 
of a chain, and N c is the average bead number between 
consecutive entanglements. The reptation theory [2-7] is 
the most successful approach to date in describing the 
dynamics of entangled polymer chains in a surprisingly 
simple manner. It has been supported by a number of 
experimental [8,9] and numerical [10-20] papers, where 
experimentally accessible quantities such as the stress 
relaxation function G(t), the incoherent dynamic scat- 
tering function, the mean square displacement, and the 
viscosity have been compared with the theoretical predic- 
tions [3]. However, N c estimated from the mean square 
displacement and that from the plateau modulus were 
around 30 and 70, respectively [13], where the difference 
of these two estimations indicates inaccuracy of the pref- 
actors in the predicted formulae [3,21]. These simula- 
tions were not direct observations of entanglements, but 
rather confirmed the existence of the entanglement ef- 
fects indirectly using the formulae of the reptation theory. 
Measurements of the dynamic structure factor S(k, t) by 
the neutron spin-echo method [8] gave information of the 
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tube diameter oc N e ' in the reptation theory and the re- 
sultant N e was consistent with the estimated value from 
the plateau modulus (both being of order 100). 

Understanding macroscopic rheological properties of 
polymer melts in terms of microscopic molecular dynam- 
ics is also of great importance, where the mechanism of 



shear-thinning behavior is very different depending on 
whether N < N c or N > N c and whether the glass tran- 
sition is approached or not. Though still inadequate, ex- 
tensive simulations have been performed in this direction 
[15-18,22-24]. 

The reptation theory [2,3] treats entanglements as dis- 
crete objects severely constraining the chain motions only 
through "tubes". However, the exact nature of entangle- 
ments is not yet known and it is highly nontrivial how 
to detect them "directly" in simulations. Therefore, it is 
at present impossible to determine N c "precisely" from 
simulations. We also point out that the reptation the- 
ory does not provide the scaling functions of the physical 
quantities applicable even for not large N/N e . The main 
aim of this paper is hence to give attempts to identify and 
visualize entanglements on the microscopic level. Here 
we mention related simulations. Ben-Nairn et al. [25] 
could visualize individual entanglements using the fact 
that incidental contacts of the particles and entangle- 
ment contacts behave differently because entanglement 
constraints are long-lived. Gao and Weiner introduced a 
time-averaged atomic mobility and found that intrachain 
atoms of relatively low mobility tend to cluster in group 
along the chain for N = 200, which suggests the existence 
of entanglements [14]. In recent Kroger and Hess' sim- 
ulation for 10 < N < 400 [17], the zero-shear viscosity 
rj (obtained at extremely weak shear) changed over from 
the behavior rj oc N to the behavior rj oc iV a with a in the 
range of 3 and 3.5 around N = iV c ~ 100. This crossover 
polymerization index A^ c should be comparable to N c but 
they did not set N c /N c = 1, probably because of the lack 
of the theoretical scaling formula rj — Nf v \ s (N/N e ) de- 
scribing the Rouse-to-reptation crossover as a function of 
N/N e . In a similar simulation, Aoyagi and Doi [18] cal- 
culated the steady state viscosity and normal stress dif- 
ferences to examine nonlinear rheology for N = 100, 200, 
and 400. 

In Table I, we summarize the simulations of freely 
jointed bead-spring (Kremer-Grest type) chains. The es- 
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timated values of N e are only crude ones. Except for 
ours in this paper, they were obtained indirectly by fit- 
ting of numerical data to the predictions [21]. At present 
it is still difficult to perform sufficiently large simulations 
with N > N c . 

The organization of this paper is as follows. In IIA 
our model system and our numerical method will be ex- 
plained. In IIB, in the range N = 10 — 250, we will exam- 
ine the time-correlation function of the end-to-end vector 
C(t), the stress relaxation function G(t), the zero-shear 
viscosity obtained from rj = J °° dtG(t), and the diffusion 
constant D. These data will indicate N c ~ 100 with the 
aid of the reptation theory in agreement with the previ- 
ous simulations. The rheological crossover from the New- 
tonian to shear-thinning behavior will also be examined 
for various N as in the work by Aoyagi and Doi [18]. We 
here stress that G(t) exhibits multi-scale relaxations over 
many decades in chain systems and its numerical calcu- 
lation has been rare because it requires very long simula- 
tions [19,20,23,24]. We will then present time-averaging 
methods of detecting entanglements in quiescent states 
in IIC and in rapidly sheared states in IID. Our "direct" 
observations of entanglements will again give N c ~ 100 
and enable us to examine how the stress overshoot and 
chain stretching occur in transient states [6,26,27]. 



II. NUMERICAL SECTION 
A. Model 

We used the bead-spring model [10] for our polymer 
melts composed of M chains with N particles or beads 
in a cubic box with volume V. The total particle number 
NM was 1000 for N = 10, 25, and 100, and was increased 
to 2500 for N = 250. All the particle pairs interact via a 
truncated Lcnnard-Jones potential defined by [10] 

U u {r)=Ae[{a/rf 2 -{a/rf]+e (1) 

The right hand side is minimum at r = 2 1 / 6 er and the po- 
tential is truncated (or zero) for larger r. Using the repul- 
sive part of the Lennard- Jones potential only in this man- 
ner, we may prevent spatial overlap of the particles [10]. 
The number density was fixed at n = NM/V = l/er 3 , 
and the temperature was kept at T = e/ks- At this 
temperature there was no glass-like enhancement of the 
structural relaxation time, but at T = 0.2e/fcs the 
present model with N = 10 became glassy in our pre- 
vious simulation [24]. Note that our density value is 
higher than the widely used value n — 0.85/cr 3 in the 
previous simulations [10,13,15-18,25] (With increasing n 
the free volume for particles decreases and hence N c in 
our case should be somewhat shorter than in the previ- 
ous simulations). The consecutive beads on each chain 
are connected by an anharmonic spring potential of the 
form, 



^feneM = -ifccRo Ml - (r/Ro) 2 } (2) 

where k c = 30e/a 2 and R = 1.5a. The bonded pairs 
in the same chain thus interact via the sum of the two 
potentials, Ut{t) = C^lj(?*) + C^fene('"), which has a deep 
minimum at r = 6 m ; n = 0.96cr. In our simulations the 
actual bond lengths between two consecutive beads re- 
mained very close to b min with deviations being at most 
a few % of b m i n even under rapid shearing. In fact, the 
expansion /7r(r) — C/T(&min) — 478e(r/6 m ; n — l) 2 follows 
around the minimum, so the deviation of this potential 
from the minimum value becomes of order e even for 
t — frmin ~ 0.04(T. (This means that the thermal fluc- 
tuation of the bond lengths is of order O.Q4a(kBT/e) 1 ' 2 .) 
Hereafter we will measure space and time in units of a 
and To = (mcr 2 /e) 1 / 2 , respectively, with m being the par- 
ticle mass, unless confusion may occur. We numerically 
solved Newton's equations of motion and took data after 
long equilibration periods of order 10 6 . 

In quiescent cases, we imposed the micro-canonical 
condition with time step At = 0.005 under the periodic 
boundary condition. In the presence of shear flow, we 
set At = 0.0025 and kept the temperature at e/fce using 
the Gaussian constraint thermostat to eliminate viscous 
heating [28,29] . After a long equilibration period in a qui- 
escent state in the time region t < 0, all the particles ac- 
quired a velocity jy in the x direction at t = 0, and then 
the Lee-Edwards boundary condition [28,29] maintained 
the simple shear flow. The periodic boundary condition 
was imposed in the x direction. Steady sheared states 
were realized after transient behavior. 

For the case N = 250 the system length V 1 / 3 is 
2500 1 / 3 = 14, which is of the same order as the end-to- 
end distance of the chains (~ WV 1 / 2 ). For such a small 
system size under the periodic boundary condition, how- 
ever, it is not clear how accurately we can simulate the 
dynamics of real entangled polymers. In future work, 
simulations are desirable in larger systems where both 
N > Ac and V 1/3 > aN 1 / 2 are well satisfied. 

B. Crossover from Rouse to Reptation Dynamics 

First we show that our numerical results near equilib- 
rium are consistent with the Rouse or reptation theory 
[3]. In Fig.l, for various A on a semi- logarithmic scale, 
we show the normalized time-correlation function of the 
end-to-end vector P = Rjy — Rq, 

C(t) = {P(t + t Q )-P(t Q ))/{\P(t )\ 2 ) (3) 

which is normalized such that (7(0) = 1. Here, because 
of the finite size of our system, the denominator and nu- 
merator on the right hand side remain to be considerably 
dependent on t even after taking the averages over all 
the chains [30]. Thus they were also averaged over the 
initial time to- The statistical (temperature-dependent) 
bond length b is defined by 
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(\P(t a )\ 2 ) = (N-l)b 2 (4) 

In the present case b is equal to 1.25a and is longer 
than the actual bond lengths \bj\ = \Rj + i(t) — Rj(t)\ = 
b m in(~ 0.96a), where 6 m i n was introduced below eq.2 [31]. 
For this time-correlation function both the Rouse dynam- 
ics and the reptation dynamics predict the same simple 
functional form [3] . For N 3> 1 it is a scaling function of 
t/r T in terms of a relaxation time r r as 

C W= E ^«P(V^r) (5) 
odd p 

where the summation is over odd p and 1 < p < N — 1. 
Since the first term (p = 1) in the summation is dominant 
for any t,@ this function decays nearly exponentially and 
r r may be determined from C(r r ) = e . As shown in 
Fig. 1, for N = 10, 25, 100, and 250, we obtained 

r r = 270, 1850,4.6 x 10 4 , 6 x 10 5 (6) 

or T r /N 2 ~ 2.7,3.0,4.6, and 9.6, respectively. The calcu- 
lated curve for N = 100 can be excellently fitted to the 
theoretical function C{t) in cq.5. For the other values of 
N the deviation is at most of order 4% around t ~ 0.1r r . 
However, at long times t > r r , good agreement between 
the calculated and theoretical curves was obtained for 
any N. 

Theoretically [3] , the relaxation time r r should be equal 
to the Rouse relaxation time 

tr = Toi TV 2 (7) 

for N N c and to the reptation or disentanglement time 
[3] 

r d = r^N 3 /N e (8) 

for N > N c . Here r i = Qb 2 /3ir 2 k B T and = 
Q) 2 /ir 2 k B T in terms of the monomeric friction constant 
£ and the bond length b in eq.4. If £ and b are assumed 
to be independent of N, the theory predicts 

r' m = 3r i, r d /r R = 3N/N C (9) 

It has been argued that the Rouse time tr has a well- 
defined physical meaning even in the reptation regime as 
the relaxation time of the chain contour in a tube [3,6,7]. 

For N = 10, the Rouse dynamics should be valid and 
T r = tr should hold, so our data of b and r r yield 

Toi = 2.7, (b 2 /k B T = 80 (10) 

For the other values of N we have tr = 0.0277V 2 . Our 
value of £ is about twice larger than in the previous 
simulations with na 3 = 0.85 (see Ref. [13]). We then 
discuss the case of N = 250, for which our data give 
ty/tr - 4.6 x 10 4 /(0.027 x 250 2 ) 3.6 and the theoret- 
ical results in eq.9 give Td/TR = 750/N c . In this paper 
(from Figs. 2, 5-7, and 10 below), we will obtain N c ~ 100, 



which then yields t^/tr ~ 7.5 and t^)t t ~ 7.5/3.6 ~ 2. 
Thus there arises a difference of a factor 2 between the 
theoretical Td and the numerically obtained r r . It could 
stem from two possible origins. One is that the Rouse- 
to-reptation crossover has not yet been well realized for 
N = 250 ~ 2.5N C . The other is that the prefactors in 
the reptation theory are inaccurate as suggested by Piitz 
et al. for large N [13]. These seem to be both relevant 
in the present work. 

In Fig. 2, we show our numerical data of the stress re- 
laxation function for N — 25 and 250 expressed by 

G(t) = {kvT)- l V{a xy (t + to)a xy (t )) (11) 

where V is the system volume and the tensor a a p(t) is 
defined by [15,16',24] 

V 1 J drU a0 (r, t) = p(t)6 afj - a a[i (t) (12) 

in terms of the microscopic stress tensor H a ^(r,t). The 
pressure p(t) may be defined such that a a p(t) is traceless 
or deviatoric. See Refs. [15,16,24] for the microscopic ex- 
pression of a a fj (t) . On one hand, the curve for the short 
chain case N = 25 and M — 40 is a result of the aver- 
ages over 10 independent runs and over the initial time 
to. For t > 10 it can be excellently fitted to the Rouse 
relaxation function (left dotted line) [3], 

Gn(t) = ^ J2 ex P ( Vt/roiW 2 ) (13) 
P =i 

where toi = 2.7 as determined in eq.10. However, for 
the longest chain case N — 250 and M = 10, G(t) was 
calculated from a single very long run performed up to 
t = 5 x 10 6 = 10 9 Ai with the average over to being 
taken. For 10 < t < 10 3 , it agrees with G R (t) (right dot- 
ted line), but in the terminal time range ( > r r = 6 x 10 5 
the calculated G(t) relaxes much slower than predicted 
by the Rouse dynamics. In Fig. 2, although the charac- 
teristic plateau behavior for N 3> N c is not clearly seen, 
we write the theoretical stress relaxation function G rep (t) 
given by [3] 

G rep (t) = ^C(t) (14) 

with iV = 100 (dashed line), where C(t) is defined in 
eq.5. The agreement appears to be fair for t > 10 5 , but 
due to the noisy behavior and the absence of well-defined 
plateau in G(t), we may claim only that N e is in the range 
70 < N c < 150 from the fitting. Note that the value of r r 
at N = 250 used for C(t) in cq.14 is obtained from Fig.l 
or eq.6. 

Further remarks regarding Fig. 2 are as follows, (i) 
First, the initial value G(0) = V {a 2 y ) /ksT takes a very 
large value about 92 (not seen in Fig. 2) and is nearly 
independent of N. Because the initial values of Gr(£) 
and G rcp (t) are much smaller as Gr(0) = nk B T <~ 1 
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and G rop (0) ^ 3nk B T/5N e ~ 0.25, agreement of G(t) 
with these model relaxation functions is attained only 
after transient relaxations [32]. As reported previously 
[19,23], the high-frequency vibration of the bond lengths 
(with period 0.14 here) gives rise to initial oscillatory be- 
havior in G(t) at short times (t < 1) [33]. (ii) Secondly, 
we notice marked noisy behavior of the curves in Fig. 2 
for t > t T as already reported in Refs. [19,24]. For such 
large time separation, the correlation of order 10 -4 of the 
initial value needs to be picked up, while the amplitude of 
the thermal fluctuations of <r xy (t) at each time t is given 

by 



'xy 



-1 1/2 



[G(0)k B T/V] 



1/2 



(15) 



This quantity is of order (92/2500) x l 2 ~ 0.2 for the case 
of N = 250 in the dimensionless units [24]. The noisy 
behavior in G(t) in the terminal time range can in prin- 
ciple be eliminated only if the system size is very large 
and/or many runs are performed. 

In Fig. 3, we summarize our results of the relax- 
ation time r r , the zero-frequency shear viscosity r\ = 
J °° dt G(t) in the linear regime, and the diffusion con- 
stant D, as functions of N. Here r) oc N in the Rouse 
dynamics and 77 oc N 3 /N 2 in the reptation dynamics. 
The diffusion constant was obtained from the relation, 



(\AR G (t)\ 2 )=6Dt 



(16) 



where t > r r and ARc(t) — Rc{to + 1) — i?c(io) is the 
displacement vector of the center of mass of the chains 
in a time interval [t ,t + to] with width t. It is known 
that D oc TV -1 in the Rouse dynamics and D oc N e /N 2 
in the reptation dynamics [3]. In Fig. 3 we can see that 
the dynamics of our system cannot be described by the 
Rouse dynamics with increasing N(> 100). Though the 
largest N(— 250) is only 2.5 times larger than our es- 
timated N e , the N dependencies in Fig. 3 are consistent 
with the crossover from the Rouse to reptation dynamics 
reported in the previous simulations [10-13,15-17]. 

We also show that our model polymer melts exhibit 
strongly nonlinear behavior for 7 > r r _1 . Fig. 4 displays 
the steady-state shear viscosity 77 in (a) and the steady- 
state normal-stress difference Ni in (b), where 



77(7) = {^x V )h, #1(7) = - (vvv) 



(17) 



Here the stress components consist of the contributions 
from all the particles as defined by eq.12 and we took 
the time-average (ct^ --(to)) over to within an interval with 
width 200. The arrows indicate the points at which 
7 = t~ . The 77(7) at high shear is nearly independent of 
N in agreement with rheology experiments [1,7,34] and 
the previous simulations [15-18], while Ni(j) is an in- 
creasing function of N at any shear rate. If the data of 
N = 100 and 250 at high shear are fitted to power laws, 
we roughly obtain 77(7) ~ A f~ a and ^1(7) ~ 7 ai iV Cl with 



a ~ 0.7, ai ~ 0.5, and c\ ~ 1.0. The second normal 
stress difference ^2(7) = (cj/j/) — (&zz) (not shown here) 
was smaller than ^1(7) by about two orders of magni- 
tude at any shear rates. In the literature [3,6,7], it has 
been argued that the relaxation of the chain contours oc- 
curs on the time scale of tr and is of crucial importance 
in nonlinear rheology under rapid deformations. In our 
case, however, Fig.l shows that the ratio tv/tr is only 
about 3.6 even for N — 250, so we cannot draw definite 
conclusions on the overall nonlinear rheology 



C. Entanglements in Quiescent States 

For the longest chain case N — 250, we attempt to 
identify and visualize entanglements. We expect that 
there should be singular enhancement in the Lennard- 
Jones potential energy between particles near an entan- 
glement point. To examine this effect, we first define the 
potential energy of non-bonded interaction on the i-th 
particle by 



W i (t)= J2 U LJ (\r t (t)-r,(t)\) 



(18) 



j'Gnon— bond 



The particle pairs i and j mostly belong to different 
chains giving rise to the inter-chain contributions in 
eq.15. However, we also include the contributions from 
the pairs belonging to the same chain but being not ad- 
jacent to each other (j ^ i ± 1). We found that the 
thermal fluctuations of Wi (t) are so large that their dis- 
tributions are nearly Gaussian at each time t without 
any noticeable correlations even between adjacent beads 
(i and i + 1 on the same chain). The non-bonded inter- 
actions Wi(t) consist mostly of rapidly varying thermal 
fluctuations uncorrclated to one another. With the end of 
reducing such rapid components, we introduce the time 
average of Wj(rj), 

Wi(t,r) = - [ dt'Wi{t + t') (19) 

T JO 

where the time interval r is much longer than the micro- 
scopic time To(= 1 in our units). This time average is 
analogous to that introduced by Thirumalai and Moun- 
tain in analyzing dynamics of supercooled liquids [35]. 
We introduce the variance cr(r) [35], 



^) 2 = J^HlWi(to,r)-(W)] 2 

i 

= 1 f dt'(l-t'/T)F NB (t') (20) 
T Jo 



where 



F (f') = (Wi(t + t')Wi(to)) - {Wf 



(21) 



Here (W) is the thermal average of Wi(to) over all the 
beads nearly independent of i and t ■ In the second line of 
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eq.20, the time-correlation function F NB (t') is assumed 
to be independent of the initial time to- We found that 
ct(t) 2 decays nearly as At^ 1 for t>1. The coefficient A 
should then be given by A = 2 dt'F NB (t'). This indi- 
cates that most of the contributions in the time integral 
in eq.18 behave as thermal white noise [35]. 

From the reptation theory [3,8] the distance of the ther- 
mal bead motions during a time interval of r is estimated 
in the short-time range r <C r c as 

£(t) = v/(|A#(t)| 2 ) ~ a^/re) 1 / 4 ~ r 1 ^ ( 22 ) 

where AR(t) = R(t n + r) — R(t ) is the displacement 

1 /2 

vector of a bead during a time interval of r, a = bN c ' 
is the tube diameter, and [36] 

t c ~ t i7V 2 - 3 x 10 4 (23) 

is the onset time of the effect of tube constraints. The re- 
lation in eq.22 is well satisfied in the range 10 < t < 10 4 
for N = 250 in our simulation. This power law has 
been confirmed in the previous simulations [10,11,13]. To 
achieve visualization of entanglements in the following, 
we should require I < a ~ A^ 2 and hence r < t ~ TV 2 
in eq.19. 

With the above time-averaging procedure, the white 
noise should be mostly eliminated with increasing r and, 
as a result, long-lived correlations due to a small num- 
ber of entanglements should become detectable in the 
range 1 r < t c . In order to demonstrate this, 
in Fig. 5, we display normalized instantaneous values 
\Wi{t) — (W)]/a(0) in (a) and normalized time-averaged 
values [Wi(t,r) - (W)]/ct(t) for r = 5 x 10 3 = 0.8 x 
10~ 2 T r ~ 0.2t o in (b) at an appropriate time t after a 
long equilibration period. For this r, the distance £(t) in 
cq.22 is given by 4.5. Here the 10 chains in our system at 
time t are straightened horizontally for the visualization 
purpose (Similar pictures of a time-averaged atomic mo- 
bility were given by Gao and Weiner [14]). In (a) almost 
no correlation can be seen along the chains. In (b), on 
the other hand, several consecutive beads form "active 
spots" having relatively large values of Wi(t,r). Here a 
bead (being the i th one in a chain) is defined to be active 
if the average 

, i+3 

W i (t,r) = - ]T W s {t,T) (24) 

j=i-3 

over the 7 adjacent values (j = i — 3, i — 2, • • • , i + 3) in 
the same chain is larger than l.ler(r). This averaging 
"in space" furthermore eliminates random, small-scale 
fluctuations consisting of a few beads and, as a result, 
the variance of the average Wi(t,r) becomes 0.65cr(r) 
for t = 5 x 10 3 . Then we select 4% of the total beads as 
active ones. In Fig. 5b they form the active spots num- 
bered from 1 to 18, consisting of several consecutive ac- 
tive beads. In Fig. 5c we show Wi(t, r) for the three chains 
with the active spots from 4 to 7. Here we are expecting 



that these active spots should arise from entanglements 
in most cases except for accidental enhancement of the 
non-bonded interactions. Fig. 5b indicates the existence 
of two or three entanglements on each chain leading to 
the estimation iV e = 100. 

To examine the correlations in Wi(t, r) along the chain 
contour quantitatively, we define the intra-chain correla- 
tion function of the non-bonded interaction 

E{n ' T) = W) ? [Wi(t,r) - (wmwjfrr) - (w)] 

j^i—n 

(25) 

where the two beads i and j are separated by n on 
the same chain and the average over all the chains is 
taken. The normalization factor 7V(t) is defined such 
that E(0, t) = 1; then, 7V(t) = MNa(r) 2 oc r. In Fig.6, 
we show E(n, t) in (a) and its Fourier transformation 

JV-l 

P(k, T ) = E ( n ' T ) cos(fcn) (26) 

71=0 

in (b) for r = 0, 50, 500, 5 x 10 3 , and 5 x 10 4 . The longest 
r is of the same order as r e in eq.23. The most conspic- 
uous feature is that E(n, r) takes a negative minimum 
around n = 45 and positive maxima around n = 80 with 
increasing r. The average displacement i{r) in eq.22 was 
calculated to be 4.5 for r = 5 x 10 3 and 7.7 for r = 5 x 10 4 . 
The contributions giving rise to these extreme grow in 
time in the normalized correlation E(n, r) or equivalently 
decay slower than r _1 in the unnormalized correlation 
E{n, r)7V(r). Correspondingly, the Fourier transforma- 
tion P(k,r) has a peak at k = 27r/90 at large r. This 
suggests N c = 90. In addition, if r is much larger than 
t c , no periodic structure was observed in E{n,r) (not 
shown here). This should be because the entanglements 
are delocalized along the chains in long time intervals 
with r 3> t c . Furthermore, as can be seen in the inset 
of Fig. 6a, the correlations between nearby beads E(n, r) 
with 1 < n < 10 are nearly zero for r = but increases 
with increasing r. We may determine the characteristic 
width n w (r) by E(n,r) > 10~ 2 for n < n w (r). If this 
definition is used, the calculated values of n w (r) and £(t) 
in eq.22 nearly coincide (within a few 10%). 

D. Entanglements under Rapid Shearing and Stress 
Overshoot 

Next, we applied a shear flow with rate 7 = 10~ 3 ~ 
600/r r - 170/t r to the system of N = 250 and M = 10. 
We used the same initial values for the particle positions 
and momenta as those which produced the data shown in 
Fig. 5b. This is convenient to examine how entanglements 
behave in the quiescent and sheared conditions starting 
at exactly the same conditions. Here the time scale of 
the flow-induced chain deformations (~ 7 _1 ) is much 
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shorter than t c in eq.23. In Fig. 7 the chain conformations 
are projected onto the xz plane (perpendicular to the 
velocity-gradient direction) at jt = 5 in (a) and jt = 10 
in (b). In (b) the shear stress takes a maximum as will 
be shown in Fig. 9. Because the chains are rapidly elon- 
gated, they eventually take zigzag shapes bent presum- 
ably at entanglements. The non-bonded interactions in 
these zig-zag points become increasingly amplified with 
increasing strain. This should be because a considerable 
fraction of the stress is supported by entanglements in 
strong deformations. As a result, active spots in Wi(t, r) 
can be detected even with much smaller r than in the qui- 
escent case, so we set r = 500 = 0.5/ j at this shear rate. 
In our simulation the noise effect in Wi(t, r) is much more 
reduced for r ~ 7 -1 than in the quiescent case. Roughly 
speaking, a 2/3 fraction of the numbered active spots 
without shear in Fig. 5b remain to be active spots under 
shear strain of 0.5, and a 1/2 fraction of them become 
bent under shear in Fig. 7a and Fig. 7b. As in the crite- 
rion in Fig. 5, the definition of the active beads is given 
by Wi(t, t) > 1.1ct(t) (but with much smaller r) and the 
number of the active beads is 4% of the total bed number. 
We assign the same numbers to these active spots if their 
contour distance between the locations along the chain in 
the quiescent and sheared cases remains shorter than 10. 
The bend regions marked by + in Fig. 7, however, do not 
correspond to the numbered hot spots in Fig. 5. 

We can also see that the number of the bends has 
not decreased from Fig. 7a to Fig. 7b, but several of them 
are approaching the chain ends and will disappear (not 
shown here) . Notice that the shear stress is maximum at 
the time of Fig. 7b. In the reptation theory it is assumed 
that entanglements can be released only when they reach 
a chain end. If our bends represent entanglements, the 
disentanglement process induced by shear flow is going 
to start in Fig. 7b, then leading to a decrease of the shear 
stress. 

For comparison, in Fig. 8 we show snapshots of the 
chains for the shorter chain case of N — 25 in a quies- 
cent state in (a) and under shear 7 = 10~ 2 in (b), where 
r r = 1850(= tr) and there is no entanglement. The non- 
bonded interactions are written with r = 50, but we can- 
not see any particular meaning in its heterogeneities on 
the chains with and without shear. In strong shear flow 
7Tr ^> 1, the chains are at most times elongated along 
the flow but undergo random tumbling motions with pe- 
riod much longer than 7 -1 . In (b) compactly shaped 
chains are in the course of tumbling [24] . 

In the case of entangled melts, the shear stress cr xy {t) 
and the normal stress N\{t) = cr xx (t) — <J yy (t) are 
known to exhibit overshoot behavior in rapid shearing 
[6,7,15,16,18,27]. In Fig. 9a we show their time evolution 
after application of shear at t = with 7 = 10~ 3 in 
the run of Fig. 7. On one hand, cr xy (t) takes a maximum 
at jt = 10, at which the disentanglement starts. On 
the other hand, N\{€) takes a maximum at jt = 20 af- 
terwards. The same tendency of successive maxima of 
Vxyit) and Ni(t) was predicted theoretically [7] and ob- 



served experimentally under rapid shearing [6,26,27]. In 
Fig. 9b we show the numerical result of the orientational 
tensor Q a p{t) defined by 

N-l 

= M £ ATT T ^ & minM*)M*) (27) 
chain j— 1 

where b^ n bj(t) are the normalized bond vectors and 
^2 a Qaa — 1 since \bj\ = b m i n = 0.96ct as stated be- 
low eq.2. For this shear rate (7 = 10~ 3 ), Q xx {t) — Q yy (t) 
remains considerably smaller than 1 and the chain bonds 
are still weakly oriented along the flow direction. The 
overshoot of Q xx (t) — Q yy (t) indicates retraction of the 
chain contours (tubes) after onset of disentanglement. 
Comparing Fig. 9a and Fig. 9b, we notice the proportion- 
ality of the two deviatoric components of the two tensors 

a xy {t) = A Q xy (t), JVi(t) = A [Q xx (t) - Q yy (t)} (28) 

where Aq = 2.2 (with the stress components being mea- 
sured in units of ecr~ 3 ) [37]. Note that the deviatoric 
part of the stress of polymer melts is believed to be 
nearly equal to that of the entropic stress contribution 
(<~ k^TnQ a (i) far above the glass transition temperature 
[2,3]. Furthermore, the deviatoric part of the dielectric 
tensor e Q) 3 is proportional to that of Q a /3 provided that 
the microscopic polarization tensor is uniaxial along the 
bond direction. Thus we obtain the well-known stress- 
optical relations 

^xy — Cq(J ' xy , £ xx ^yy Co((T xx ^yy) (29) 

where Co is a polymer-dependent constant. In Fig. 9c we 
also show the following angle 

X = \tnn- l {2<j xy /N 1 ) (30) 

On the basis of the stress-optical law, this angle is mea- 
sured as the extinction angle \ = (1/2) tan -1 [2e xy /(e xx — 
e yy )] in birefringence experiments. In Fig. 8c we can see 
that % exhibits a small undershoot around jt = 30 after 
the peaks of o xy (t) and N±(t). A similar retarded un- 
dershoot was observed experimentally but has not been 
explained theoretically [7]. 

As the final example, we examine the case of much 
larger shear rate 7 = 10~ 2 - 6000/r r - 1700/r R . Fig.10 
displays the snapshots of the chain conformations for 
jt = 5 in (a) and for jt = 12.5 in (b), where the 
chain stretching is stronger than in Fig. 7 and the num- 
bers of entanglements remain unchanged. The time in- 
terval r in eq.16 is set equal to 50. In Fig. 11 we show 
&xy{t) and N\(t) in (a), both exhibiting a peak around 
jt = 12.5, and Q xy (t) and Q xx (t) — Qyy{t) in (b). The 
stress overshoot is more enhanced than in the smaller 
shear case in Fig. 9, and the stress components decrease 
rather abruptly with onset of disentanglement. However, 
Qxx{t) — Q yy {t) saturates to a value about 0.7 without ex- 
hibiting overshoot. Here even the bonds themselves align 
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in the flow direction and the chain stretching becomes 
nearly complete. Interestingly, this bond alignment is 
still maintained even after onset of disentanglement. As 
we remarked below eq.2, bond elongation of order 3% is 
anharmonic for the potentials in eq.l and eq.2 and gives 
rise to a tensile force of order e/cr. In this simulation, 
such strong forces are exerted on most of the bonds and 
the proportionality relation in eq.25 does not hold. 



those in the experiments [6,26,27]. 
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III. SUMMARY AND CONCLUDING REMARKS 



We have presented attempts to detect and visualize the 
entanglements in a model polymer melt together with at- 
tempts of indirectly deriving N c as in the previous sim- 
ulations. All the methods have yielded N e ~ 100. We 
admit that the visualization in the quiescent case is not 
yet firmly established by themselves in view of the fact 
that the thermal noise still affects the data even after 
averaging in space and time as in Fig. 5c. However, un- 
der rapid shearing, a large fraction of active spots with 
relatively large nonbonded interactions become bent, ev- 
idently indicating the existence of obstacles for the chain 
motion. Remarkably, the active spots in the quiescent 
and sheared cases in Fig. 5 and Fig. 7 coincide with a large 
probability (~ 2/3). We claim that a large fraction of 
such obstacles arise from entanglements preexisting even 
before application of shear. However, a few bends and 
hot spots in Fig. 7 are not detected by the visualization 
method in Fig. 5. This suggests that some hot spots from 
our method may not represent entanglements. 

As remarked already at the end of IIA, the system size 
in our simulations is still comparable to the end-to-end 
distance for N = 250 and our results need to be fur- 
ther checked in future large-scale simulations with longer 
chains. To eliminate the large thermal noise effect in the 
quiescent case, we should take data in the well-defined 
reptation regime under N ^> N c . 

Performing very long simulations, we have also cal- 
culated the stress relaxation function G(t), which ex- 
hibits the Rouse-to-reptation crossover with increasing 
the polymerization index N, and studied nonlinear rhe- 
ology in transient and steady sheared states. In tran- 
sient states under shear, the stress overshoot sets in as 
the bends approach the chain ends and disappear, as can 
be seen from Fig. 7b and Fig. 9. This is also one of our 
main results giving molecular information on the stress 
overshoot under rapid shearing. 

In real long chain systems, the ratio t^/tr ~ N/N e can 
be very large. Hence, in shear flow, there can be three 
characteristic shear regions [7] given by (i) 7 < tT 1 , (ii) 
< 7 < Tj^ , and (iii) 7 > . Nonlinear shear ef- 
fects emerge in the regions (ii) and (iii), while the linear 
response theory in terms of G(t) in cq. 1 1 is valid only in 
the region (i). In our study, the intermediate region (ii) 
is not well-defined, but the calculated overshoot and un- 
dershoot relaxations in Fig. 9 (in the region (hi)) resemble 
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TABLE I. Summary of simulations of freely jointed chains. The second column shows the chain length in the simulations. 
The third and fourth columns show the system volume V and the simulation time £ max in units of the chain dimension bN 1 ^ 2 
and the Rouse time tqN 2 . In the fifth to tenth columns the estimated values of N e are given (if estimated). 
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FIG. 1. The normalized time-correlation function C{i) for 
the end-to-end vectors in eq.3 for various (solid lines) on a 
semi-logarithmic scale. They can be fitted to the theoretical 
expression in eq.5 (dashed lines). The arrows indicate the ter- 
minal relaxation time r r defined by C(r r ) = e _1 , which is the 
Rouse relaxation time tr or the reptation time t<j depending 
on whether N < N a or N > N e . 
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FIG. 2. The stress relaxation function G(t) for N — 25 
and 250 (solid lines) in units of ecr~ 3 . Its N dependence is 
weak for t < 10 2 but becomes increasingly stronger at later 
times. For N — 25 the curve of G(t) can well be fitted to the 
Rouse relaxation function Gr(£) with roi = 2.7 and N = 25 
in eq.10 (left dotted line). For N = 250, it can be fitted to 
the reptation relaxation function G Tep (t) in eq.12 for t > r r 
(dashed line), where r r = 6 x 10 5 in C(t) and A^e = 100 in the 
prefactor. For comparison we also show Gn(t) with roi = 2.7 
and AT = 250 (right dotted line). It much deviates from the 
calculated curve of G(t) of N = 250 for t> 10 4 , while they 
agree fairly well for 10 < t < 10 3 . 
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FIG. 3. iV-dependence of N 2 /t t (•), N/r) (■), and ND 
(O). These quantities are expected to tend to the Rouse 
limits for 1 < iV < ft. The solid lines represent 
(3n 2 k B T/(b 2 )/(l +3N/N e ), (36/n(b 2 )/(l + 8N 2 /5N 2 ), and 
{k B T/C)/(l + 3N/N e ) with iV c = 100. These approximate ex- 
pressions extrapolate the predicted formulae 3 of the Rouse 
and reptation behaviors in the two limits N <C N c and 
N > N c . 
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(b) Non-bonded interactions with t = 5000 




FIG. 5. Distributions of the non-bonded interaction energy on the chains for N = 250 in a quiescent state. The 10 chains 
are straightened on the plane. In (a) the normalized values [Wi(t) — (W)]/a(0) are shown, where r = and no correlations 
along the chains can be seen. In (b) the normalized, time-averaged values \Wi{t,r) — (W)]/a(r) with r = 5x 10 3 are shown, 
which are distinctly large in line segments consisting of several beads (in orange) presumably due to entanglements. In (b) the 
active spots are numbered from 1 to 18 according to the criterion given around eq.24. In (c) the data of [Wi(t,r) — {W)]/a(r) 
for the three chains with the spots 4-7 in (b) are shown. The horizontal axis denotes the bead numbers 1 < N < 250 for the 
three chains. The beads above the broken line are defined as "active" beads. 
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FIG. 6. The intra-chain correlation function E(n,r) of the non-bonded interactions defined by eq.25 in (a), and its Fourier 
transformation P(fc,r) defined by eq.26 in (b), for N = 250. Here r = (thin-solid line), 5 x 10 (dotted line), 5 x 10 2 (dashed 
line), 5 x 10 3 (bold line), and 5 x 10 4 (dot-dash line). With increasing r we can see development of the minimum and the 
maxima presumably due to entanglements. The inset in (a) shows E(n, r) for small n for these r. 
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N=250, y = 0.001 




i 

FIG. 7. Snapshots of the deformed chains with TV = 250 on the xz plane at jt = 5 in (a) and at jt = 10 in (b), where a shear 
flow with 7 = 10 -3 was applied at t = with the same initial chain configuration as in Fig. 5b. The non-bonded interactions 
with r = 500 are written on the chains. The active spots satisfying the criterion given around eq.24 are detected, among which 
the numbered segments correspond to those in Fig. 5b. However, the active spots which do not correspond to those in Fig. 5b 
are marked by +. The flow is in the horizontal (x) direction, and the shear gradient is in the out-of-plane (y) direction. In (c) 
the data of [VFi(£, r) — (W)]/a(r) for the three chains with the spots 4-7 in (b) are shown. 
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N=25, y = 0.01 
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(a) jt = 

FIG. 8. Snapshots of short chains with TV = 25 in the xz plane at A /t 
with 7 = 10~ 2 was applied at t = 0. 
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FIG. 9. The stress growth functions after application of 
shear flow with j = for N = 250 in (a). The arrows (a 

and b) indicate the points of jt = 5 in Fig. 7a and jt = 10 
in Fig. 7b. The maximum of <J X y{t) is at jt = 10 and that of 
Ni(t) is at jt — 20. The corresponding components of the 
orientational tensor Q a p{t) in eq.27 are shown in (b). The 
extinction angle % m eq.30 is shown in (c), which undershoots 
before reaching a larger steady-state value. These figures are 
the results of the single run displayed in Fig. 7. 
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(a) yf = 5 (b) yf = 12.5 

FIG. 10. Snapshots of the deformed chains with TV = 250 in the xz plane at jt = 5 in (a) and at jt ~ 12.5 in (b), where a 
shear flow with 7 = 10 -2 was applied at t = 0. The non-bonded interactions with r = 50 are written on the chains using the 
color map on the left. In (b) the stretching is nearly complete between the bends. 
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yt yt 
FIG. 11. The stress growth functions after application of shear flow with 7 = 1CP 2 for N = 250 in (a). The arrows (a 
and b) indicate the points of -yt = 5 in Fig. 10a and -yt = 12.5 in Fig. 10b. The corresponding components of the orientational 
tensor Q a /3{t) are shown in (b), which demonstrate bond elongation along the flow and behave very differently from the stress 
components. 
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